#include "CylinderReceiverRectangleGridRayTracing.cuh"
#include <algorithm/rayTracing/dataLoader/CloudLoader.h>
#include <algorithm/rayTracing/rayTracingDataStructure/receiver/CudaCylinderReceiver.cuh>
#include <algorithm/rayTracing/rayTracingPipeline/quasiMonteCarloRayTracing/Grid/rectGridDDA.cuh>
#include <algorithm/rayTracing/rayTracingPipeline/quasiMonteCarloRayTracing/Receiver/cylinderReceiverIntersection.cuh>
#include <algorithm/shared/tool/randomNumberGenerator/RandomGenerator.cuh>

#include <cmath>

using namespace solar::cuda;

inline __host__ __device__ auto eta_aAlpha(const float& d) -> float
{
    if (d <= 1000.0f)
        return 0.99331f - 0.0001176f * d + 1.97f * (1e-8f) * d * d;
    return expf(-0.0001106f * d);
}
/**
 * Eray = (ID ⋅ cosφ ⋅ ρ ⋅ Shsub / (Nc ⋅ Srsub) ) ⋅ ηaα
 * @param distance
 * @param dir
 * @param normal
 * @param factor
 * @return
 */

inline __host__ __device__ auto calEnergy(float distance, float3 dir, float3 normal,
                                          float factor) -> float
{
    //      cos(dir, normal)       *  ηaα                 *  (DNI * Ssub * reflective_rate /
    //      numberOfLightsPerGroup)
    return std::abs(dot(dir, normal)) * eta_aAlpha(distance) * factor;
}

__host__ __device__ static auto testdevicecuda_in_raytracing() -> bool
{
    printf("testdevicecuda_in_raytracing\n");
    return true;
}

// =========================这里的代码是从rectGridDDA移植过来的 -
// begin===========================================

template <typename T>
inline __host__ __device__ auto absDivide(const T& denominator, const T& numerator) -> T
{
    if (numerator <= Epsilon && numerator >= -Epsilon)
    {
        return T(INT_MAX);
    }
    return std::abs(denominator / numerator);
}

inline __host__ __device__ auto calTMax(float dir, float interval, int current_index,
                                        float current_pos) -> float
{
    return std::abs(float(current_index + (dir >= 0)) * interval - current_pos);
}

inline __host__ __device__ auto equal_(float n1, float n2) -> bool
{
    return (n1 < n2 + Epsilon) && (n1 > n2 - Epsilon);
}

__host__ __device__ bool less_(float n1, float n2)
{
    return n1 < n2 - Epsilon;
    //    return true;
}

//--------------------把rectGridDDA中的代码移植过来了。。。----------------------------------
__host__ __device__ static auto intersect(const float3& origin, const float3& dir,
                                          const float3* d_heliostat_vertexes,
                                          const int* d_grid_heliostat_match, int start_id,
                                          int end_id, int subheliostat_id,
                                          int numberOfSubHeliostat) -> bool
{
    // printf("===== enter intersect\n");
    float t;
    float u;
    float v;
    for (int i = start_id; i < end_id; ++i)
    {
        int subHeliostatIndex = d_grid_heliostat_match[i];
        if (subHeliostatIndex < subheliostat_id ||
            subHeliostatIndex >= subheliostat_id + numberOfSubHeliostat)
        {
            if (global_func::rayParallelogramIntersect(
                    origin, dir, d_heliostat_vertexes[3 * subHeliostatIndex + 1],
                    d_heliostat_vertexes[3 * subHeliostatIndex],
                    d_heliostat_vertexes[3 * subHeliostatIndex + 2], t, u, v))
            {
                return true;
            }
        }
    }
    return false;
}
//------------------------把rectGridDDA中的代码移植过来了。。。----------------------------
__host__ __device__ static auto collision_in_rectgrid(const float3& origin, const float3& dir,
                                                      const CudaRectangleGrid& rectGrid,
                                                      const float3* d_subheliostat_vertexes,
                                                      const HeliostatArgument& heliostatArgument,
                                                      int id) -> bool
{
    // printf("=== enter collision_in_rectgrid()\n");

    //  Step 1.1 Initial current position of origin in the scene.
    int3 pos = make_int3((origin - rectGrid.getPosition()) / rectGrid.getInterval());

    //  Step 1.2 StepX, StepY, StepZ
    int3 Step;
    Step.x = (dir.x >= 0) ? 1 : -1;
    Step.y = (dir.y >= 0) ? 1 : -1;
    Step.z = (dir.z >= 0) ? 1 : -1;

    //  Step 1.3 Initial tmaxX, tmaxY, tmaxZ
    float3 t_max; // avoid divide 0

    // tMax = |pos ± interval - (origin-rectGrid.getPosition())| / dir
    t_max.x = absDivide(
        calTMax(dir.x, rectGrid.getInterval().x, pos.x, origin.x - rectGrid.getPosition().x),
        dir.x);
    t_max.y = absDivide(
        calTMax(dir.y, rectGrid.getInterval().y, pos.y, origin.y - rectGrid.getPosition().y),
        dir.y);
    t_max.z = absDivide(
        calTMax(dir.z, rectGrid.getInterval().z, pos.z, origin.z - rectGrid.getPosition().z),
        dir.z);

    //  Step 1.4 Initial tDeltaX, tDeltaY, tDeltaZ
    float3 t_delta;
    t_delta.x = absDivide(rectGrid.getInterval().x, dir.x);
    t_delta.y = absDivide(rectGrid.getInterval().y, dir.y);
    t_delta.z = absDivide(rectGrid.getInterval().z, dir.z);

    /**
     * Step 2 - Intersection
     **/

    int3 grid_index = pos;
    int grid_address = 0;

    while (grid_index.x >= 0 && grid_index.x < rectGrid.getSubGridNumber().x && grid_index.y >= 0 &&
           grid_index.y < rectGrid.getSubGridNumber().y && grid_index.z >= 0 &&
           grid_index.z < rectGrid.getSubGridNumber().z)
    {

        grid_address = global_func::unroll_index(grid_index, rectGrid.getSubGridNumber());

        if (intersect(origin, dir, d_subheliostat_vertexes, rectGrid.getDeviceGridHeliostatMatch(),
                      rectGrid.getDeviceGridHeliostatIndex()[grid_address],
                      rectGrid.getDeviceGridHeliostatIndex()[grid_address + 1],
                      heliostatArgument.subHeliostat_id[id], heliostatArgument.helionums))
        {
            return true;
        }

        // next cell location
        bool xMinFlag, yMinFlag, zMinFlag;
        xMinFlag = yMinFlag = zMinFlag = false;

        if (less_(t_max.x, t_max.y))
        { // x < y
            if (less_(t_max.x, t_max.z) || equal_(t_max.x, t_max.z))
            { // x <= z,  +x minimal
                grid_index.x += Step.x;
                xMinFlag = true;
            }
            if (less_(t_max.z, t_max.x) || equal_(t_max.x, t_max.z))
            { // z <= x,  +z minimal
                grid_index.z += Step.z;
                zMinFlag = true;
            }
        }
        else if (less_(t_max.y, t_max.x))
        { // y < x
            if (less_(t_max.y, t_max.z) || equal_(t_max.y, t_max.z))
            { // y <= z,  +y minimal
                grid_index.y += Step.y;
                yMinFlag = true;
            }
            if (less_(t_max.z, t_max.y) || equal_(t_max.y, t_max.z))
            { // z <= y,  +z minimal
                grid_index.z += Step.z;
                zMinFlag = true;
            }
        }
        else
        { // x == y
            if (less_(t_max.x, t_max.z) || equal_(t_max.x, t_max.z))
            { // x(y) <= z, +x&y minimal
                grid_index.x += Step.x;
                grid_index.y += Step.y;
                xMinFlag = yMinFlag = true;
            }
            if (less_(t_max.z, t_max.x) || equal_(t_max.x, t_max.z))
            { // z <= x(y),   +z minimal
                grid_index.z += Step.z;
                zMinFlag = true;
            }
        }
        t_max.x += (xMinFlag) ? t_delta.x : 0;
        t_max.y += (yMinFlag) ? t_delta.y : 0;
        t_max.z += (zMinFlag) ? t_delta.z : 0;
        break;
    }

    return false;
}

// ==========================这里的代码是从rectGridDDA移植过来的 -
// end=====================================

__global__ void map_raytracing(SunrayArgument sunrayArgument, CudaCylinderReceiver cylinderReceiver,
                               CudaRectangleGrid rectangleGrid, HeliostatArgument heliostatArgument,
                               float3* d_subHeliostat_vertexe, float factor, bool res,
                               solar::CloudLoader cloud)
{

    // printf("======= in kernel func: map_raytracing  nBlocks = (%d,%d), nThreads = (%d,%d)\n",
    // gridDim.x, gridDim.y, blockDim.x, blockDim.y );
    // printf("in CUDA: factor = %f\n",factor);

    int my_id = global_func::getThreadID();
    long long number_of_micro_heliostats =
        (long long)heliostatArgument.micronums *
        (long long)heliostatArgument.helionums; // 场景内所有的微定日镜个数
    // printf("1111\n");
    // printf("myId = %d, numberOfMicroHeliostats = %d, sunrayArgument.numberOfLightsPerGroup =
    // %d\n",
    //         myId, numberOfMicroHeliostats, sunrayArgument.numberOfLightsPerGroup);
    if ((long long)my_id >=
        ((long long)number_of_micro_heliostats) * (long long)sunrayArgument.numberOfLightsPerGroup)
    { // 一个定日镜上微定日镜个数 * 一个光锥内光线个数
        // printf("enter if -- exit  %d %lld
        // %d\n",myId,numberOfMicroHeliostats,sunrayArgument.numberOfLightsPerGroup);
        return;
    }
    // printf("2222\n");
    int microsize = heliostatArgument.micronums; // 一个子定日镜上微定日镜个数
    int no =
        my_id / microsize /
        sunrayArgument
            .numberOfLightsPerGroup; // !!!!!!!!相当于每个光线都同时进行计算？？？？这就是所谓光线级别的并行吗？？？

    // int is_total_reflect=0;
    float depth = 0.003f;
    float refract = 1.523f; // 玻璃介质的折射率-----与平常说的反射率、折射率不一样
    int helio_no = my_id / sunrayArgument.numberOfLightsPerGroup /
                   heliostatArgument.micronums; // 第几个子定日镜      // ！！！！！这里不太理解
    int micro_no = my_id / sunrayArgument.numberOfLightsPerGroup %
                   heliostatArgument.micronums; // 第几个微定日镜
    int substart =
        helio_no * 3; // 为什么乘以3？？？----- 因为d_subHeliostat_vertexe中每个子定日镜只存了三个点

    // printf("helio_no=%d, micro_no=%d\n", helio_no, micro_no);
    // printf("3333\n");
    // 根据三个点计算一个子定日镜的中心和法向
    float3 vera = d_subHeliostat_vertexe[substart];
    float3 verb = d_subHeliostat_vertexe[substart + 1];
    float3 verc = d_subHeliostat_vertexe[substart + 2];
    // printf("helio_no = %d, micro_no = %d, vera %f %f %f  verb %f %f %f  verc %f %f %f\n",
    //         helio_no, micro_no, vera.x, vera.y, vera.z,  verb.x, verb.y, verb.z, verc.x, verc.y,
    //         verc.z    );

    float3 helio_pos = (d_subHeliostat_vertexe[substart] + d_subHeliostat_vertexe[substart + 2]) /
                       2; // rect对角线的中点位置
    float3 helio_norm = normalize(
        cross(d_subHeliostat_vertexe[substart] - d_subHeliostat_vertexe[substart + 1],
              d_subHeliostat_vertexe[substart + 2] - d_subHeliostat_vertexe[substart + 1]));
    if (helio_norm.y < 0)
    { // 如果法向的y小于0，表示方向反了
        helio_norm.x = -helio_norm.x;
        helio_norm.y = -helio_norm.y;
        helio_norm.z = -helio_norm.z;
    }
    // printf("helio_pos %f %f %f   helio_norm %f %f %f\n", helio_pos.x, helio_pos.y, helio_pos.z,
    // helio_norm.x, helio_norm.y, helio_norm.z);

    // printf("4444\n");
    int col_no = micro_no / heliostatArgument.col; // 该微定日镜在子定日镜中，位于第几列、第几行
    int row_no = micro_no % heliostatArgument.col;
    float3 micropos = helio_pos;   // 微定日镜的中心位置
    float3 micronorm = helio_norm; // 微定日镜的法向---先设置为子定日镜的法向
    // printf("helio_no = %d, micro_no = %d, micropos %f %f %f  micronorm %f %f %f\n",helio_no,
    // micro_no, micropos.x, micropos.y, micropos.z, micronorm.x, micronorm.y, micronorm.z);

    float error = sunrayArgument.d_samplelights[my_id % sunrayArgument.pool_size].y;
    micropos.z =
        col_no * heliostatArgument.pixel + heliostatArgument.pixel / 2 -
        heliostatArgument.sizex / 2; // 这里计算得到的micropos是在子定日镜的局部坐标系下的坐标
    micropos.y = heliostatArgument.sizey / 2;
    micropos.x = row_no * heliostatArgument.pixel + heliostatArgument.pixel / 2 -
                 heliostatArgument.sizez / 2;

    // printf("5555\n");
    //    int row_no=micro_no/heliostatArgument.col;
    //    int col_no=micro_no%heliostatArgument.col;
    //
    //    float3 micropos=helio_pos;
    //    float3 micronorm=helio_norm;
    //
    //    micropos.x = col_no * heliostatArgument.pixel +  heliostatArgument.pixel  / 2 -
    //    heliostatArgument.sizex/ 2; micropos.y = heliostatArgument.sizey/ 2; micropos.z = row_no *
    //    heliostatArgument.pixel    + heliostatArgument.pixel   / 2 - heliostatArgument.sizez / 2;
    //    printf("micropos %f %f %f\n",micropos.x,micropos.y,micropos.z);

    float3 local = micropos;
    local = global_func::local2world_rotate(local, micronorm);
    local = global_func::translate(local, helio_pos);
    micropos = local; // 把微定日镜的中心位置从局部坐标系转换为世界坐标系
    // printf("micropos %f %f %f\n",micropos.x,micropos.y,micropos.z);
    //    printf("helio_no = %d, micro_no = %d, local %f %f %f  normal %f %f %f\n",helio_no,
    //    micro_no, local.x,local.y,local.z,micronorm.x,micronorm.y,micronorm.z);
    int address = my_id % sunrayArgument.pool_size;

    // 采样数组的采样值-------采样入射光方向，并把入射光方向改为全局坐标
    //  printf("address = %d, sunray_direction = (%f,%f,%f)\n",
    //          address,
    //          sunrayArgument.sunray_direction.x,sunrayArgument.sunray_direction.y,sunrayArgument.sunray_direction.z
    //          );
    float3 dir = global_func::local2world_rotate(
        sunrayArgument.d_samplelights[address],
        -sunrayArgument.sunray_direction); // ？？？？？这里为什么要按照光线方向进行旋转？？？？
    // 入射光方向（世界坐标系）
    //    dir = -sunrayArgument.sunray_direction;
    // 计算第helio_no个定日镜的第micro_no个微定日镜的法向和坐标
    dir = normalize(dir); // 入射角dir.y>0入社向量的逆方向
    // printf("initial dir = %f %f %f\n", dir.x, dir.y, dir.z);
    float3 origin = micropos; // 微定日镜中心坐标

    // printf("in CUDA_2: factor = %f\n",factor);

    // printf("6666\n");

    if (rectGridDDA::collision_(origin, dir, rectangleGrid, d_subHeliostat_vertexe,
                                heliostatArgument, no))
    { // ！！！！！！！！！！！！！！！！！！！！！！！！！！
        // printf("exit cuda --- collision-1\n");
        return;
    }
    // printf("7777\n");
    // rectGridDDA::testdevicecuda_in_rectGridDDA();
    // if(collision_in_rectgrid(origin, dir, rectangleGrid,
    // d_subHeliostat_vertexe,heliostatArgument,no ))
    // {
    //     return ;
    // }

    // printf("8888\n");
    // bool ret = rectGridDDA::testdevicecuda_in_rectGridDDA();
    // printf("ret = %d\n",ret);

    // bool ret_1 = testdevicecuda_in_raytracing();
    // printf("ret_1 = %d\n",ret_1);

    // bool ret_2 = rectGridDDA::transfer();
    // if(rectGridDDA::testdevicecuda())
    // {
    //     printf("right!!!\n");
    // }

    float3 normal =
        micronorm; // 微定日镜法向方向，(局部坐标系，在这里还跟该微定日镜所属的子定日镜的法向一致)
    // normal=global_func::local2world_rotate(turblance,normal);
    //  采样微定日镜的法向----并转换为全局坐标
    float3 normal1 = global_func::local2world_rotate(sunrayArgument.d_perturbations[address],
                                                     normal); // 微定日镜法向方向，世界坐标系
    normal1 = normalize(normal1); // normal1是加了扰动的全局坐标系中的微定日镜的法向
                                  //    normal1=normalize(normal);

    //    printf("myId %d helioNo %d microNo %d rowNo %d colNo %d micropos %f %f %f dir %f %f %f
    //    normal1 %f %f %f\n",
    //           myId,helio_no,micro_no,row_no,col_no,micropos.x,micropos.y,micropos.z,
    //           dir.x,dir.y,dir.z,normal1.x,normal1.y,normal1.z);

    //     //先和其他定日镜求交,看是否被遮挡
    //     //如果没有被遮挡,求光线与云的交点,如果在数组中,取数组中的值,不在,取1.返回值是float值
    float3 dir_n = make_float3(0.0, 1.0, 0.0);
    float cos_theta = dot(dir, dir_n); // 这里表示入射光的逆方向与垂直方向的夹角余弦
    //    float3 dirs=dir*len;
    float3 dirs = dir; // 入射光方向的逆方向

    //=====先把云层的部分注释掉====================

    //     int cloud_type=1;//0 is online, 1 is vertical, 2 is horizontal

    //     float2 cloudstart;
    //     float len=0;
    //     if(cloud_type==0){
    //         len=cloud.start_point.y/cos_theta;  // ??????????这里表示什么？？？？？
    //         dirs=dirs*len;
    //         float times=cloud.start_point.y/-sunrayArgument.sunray_direction.y;
    //         float2 cloudcenter=
    //         make_float2(-sunrayArgument.sunray_direction.x*times,-sunrayArgument.sunray_direction.z*times);
    //         // ？？？？？这里是怎么计算的
    //         cloudstart=make_float2(cloudcenter.x+cloud.start_point.x,cloudcenter.y+cloud.start_point.z);
    //         // ？？？？？这里是怎么计算的
    // //        printf("%f %f \n",cloudstart.x,cloudstart.y);
    //     }
    //     else if(cloud_type==1){
    // //        float height=0;
    //         len=cloud.start_point.y/cos_theta;          //
    //         ------注意：定日镜的高度都设置为0，所以这里应该是，定日镜沿入射光方向到云层的距离---似乎没有考虑微定日镜与所属定日镜的高度偏差
    //         dirs=dirs*len;                              //
    //         从定日镜一直延伸到云层的向量（本质就是dir拉长了而已）
    // //        float times=height/-sunrayArgument.sunray_direction.y;
    // //        float2 cloudcenter=
    // make_float2(-sunrayArgument.sunray_direction.x*times,-sunrayArgument.sunray_direction.z*times);
    //         float2 cloudcenter= make_float2(0,0);
    //         cloudstart=make_float2(cloudcenter.x+cloud.start_point.x,cloudcenter.y+cloud.start_point.z);
    // //        printf("%f %f \n",cloudstart.x,cloudstart.y);

    //     }

    //     float2 M= make_float2(origin.x+dirs.x,origin.z+dirs.z);     //
    //     点+向量----表示一个新的点，即云层二维平面上，入射光与云层的交点

    //     float my=M.x-cloudstart.x;//M.x, M.y 是M点在世界坐标系的（x,z)坐标值    //
    //     把交点M转化成云层中的局部坐标（云层坐标系的原点为cloudstart） float
    //     mx=M.y-cloudstart.y;//需要转换成 cloud map局部坐标 int Ax=
    //     floor(mx/cloud.pixel);//A在cloud map中第几行                // 交点M在云层grid中的行列
    //     int Ay= floor(my/cloud.pixel);//A在cloud map的第几列
    // //    printf("Ax %d Ay %d  cloudstart %f %f len %f\n",Ax,Ay,cloudstart.x,cloudstart.y,len);

    //     float decay=1.0;
    //     if(Ax<0||Ax>cloud.x_size||Ay<0||Ay>cloud.y_size){       //
    //     云层划分为x_size*y_size个网格，要判断交点是否落在网格内
    //         decay=1.0;
    // //      printf("times %f cloudstart %f %f %f %f \n",times,cloudstart.x,cloudstart.y,mx,my);
    // //        printf("Ax %f Ay %f\n",Ax,Ay);

    //     }
    //     else{
    // //        printf("my %f mx %f\n",my,mx);
    // //         /*世界坐标  --> X
    // //          *         |
    // //          *         |
    // //          *         Z
    // //          * cloud map 局部坐标  --> y
    // //          *                    |
    // //          *                    |
    // //          *                    x
    // //          * M.x-cloud.start_point.x M点在局部坐标系的第几列
    // //          * M.y-cloud.start_point.z M点在局部坐标系的第几行
    // //          * A    B
    // //          *    M
    // //          * D
    // //          * 根据A,B,C,D四点的值双线性插值得到M的值*/
    // //
    //         //
    //         知道M在云层的一个网格里，但是网格顶点处才是具体的衰减系数所在的位置，要求M处的衰减系数，要对M所属的格子的四个顶点进行双线性插值

    //         int A_pos=Ay+Ax*cloud.x_size;   // 交点M所在的网格在d_cloudMap中的下标

    //         // Q代表衰减系数
    //         float Q11=cloud.d_cloudMap[A_pos];                  //A
    //         float Q12=cloud.d_cloudMap[A_pos+1];                //B
    //         float Q21=cloud.d_cloudMap[A_pos+cloud.x_size];     //D
    //         float Q22=cloud.d_cloudMap[A_pos+cloud.x_size+1];   //C
    //         //A x1,y1
    //         //B x1,y2
    //         //C x2,y2
    //         //D x2,y1
    //         float x1=Ax*cloud.pixel;
    //         float x2=(Ax+1)*cloud.pixel;
    //         float y1=Ay*cloud.pixel;
    //         float y2=(Ay+1)*cloud.pixel;

    //         decay=(Q11*(x2-mx)*(y2-my)+Q21*(mx-x1)*(y2-my)+Q12*(x2-mx)*(my-y1)+Q22*(mx-x1)*(my-y1))/cloud.pixel/cloud.pixel;
    // //        decay=1.0-decay;
    //         decay=1.0 ;
    // //        printf("%f\n",decay);
    //     }

    ////    res=1;

    if (res)
    {

        float theta = NAN; // 入射光与法向之间的夹角-----入射角
        if (std::abs(dot(dir, normal1) - 1.0F) < Epsilon)
        {
            theta = 0.0F; // 入射光与微定日镜法向重合的情况
        }
        else
        {
            theta = std::acos(dot(dir, normal1));
        }

        // 夹角小于全内反射的临界角，折射,调研得到玻璃的全内反射临界角大约为PI/5

        float3 dir_1 = -dir; // 入射方向

        float eta = 1.0F / refract; // refract_是玻璃的介质折射系数
        float k = 1.0 - eta * eta * (1.0F - dot(normal1, dir_1) * dot(normal1, dir_1));

        float theta_2 = asin(sin(theta) / refract); // 折射角
        float3 dir2 = normalize(eta * dir_1 -
                                (eta * dot(normal1, dir_1) + sqrt(k)) *
                                    normal1); // !!!!!貌似是光线进入玻璃后的折射方向----再计算一下

        float3 dir22 =
            dir2 * depth / cos(theta_2); // 折射方向，该向量就是从入射点指向光线与银层的交点
        int ran = int(sunrayArgument.pool_size / 3);
        // printf("sunrayArgument.numberOfLightsPerGroup=%d\n",sunrayArgument.numberOfLightsPerGroup);
        // 计算第一次折射向量的反射
        float3 normal2 = normalize(global_func::local2world_rotate(
            sunrayArgument.d_perturbations[(address + ran) % sunrayArgument.pool_size],
            normal)); // 银层上的法向
        // normal2=normalize(normal);
        // printf("normal %f %f %f normal1 %f %f %f normal2 %f %f
        // %f\n",normal.x,normal.y,normal.z,normal1.x,normal1.y,normal1.z,normal2.x,normal2.y,normal2.z);
        float3 dir3 = normalize(reflect(dir2, normal2)); // 银层上的反射向量
        float3 dir333 = normalize(reflect(
            dir2, normalize(normal))); // no slope     normal是没有法向扰动的法向==定日镜的主法向

        float3 dir33 = dir3 * depth / cos(theta_2); // 从银层上的点出发，指向玻璃上表面交点的向量
        float3 dir3322 = dir333 * depth / cos(theta_2); // 没有考虑slope，类似于dir33
        // printf("dir33 %f %f %f dir3333 %f %f
        // %f\n",dir33.x,dir33.y,dir33.z,dir3333.x,dir3333.y,dir3333.z);

        // 折射向量的反射vector
        float3 origin1 = origin + dir22 + dir33; // 得到从玻璃出来的点
                                                 //        float3 origin13 = origin+dir22+dir3322;

        // printf("origin1 %f %f %f origin11 %f %f
        // %f\n",origin1.x,origin1.y,origin1.z,origin11.x,origin11.y,origin11.z);

        // 实际反射的原点
        float3 normal3 = normalize(global_func::local2world_rotate(
            sunrayArgument.d_perturbations[(address + ran * 2) % sunrayArgument.pool_size],
            normal)); // 从玻璃出来的哪一点处的法向
        // normal3=normalize(normal);

        float theta_refract; // 从银层到玻璃表面的入射角（位于玻璃介质中）
        if (std::abs(dot(dir3, normal3) - 1.0F) < Epsilon)
        {
            theta_refract = 0.0F;
        }
        else
        {
            theta_refract = std::acos(dot(dir3, normal3));
        }

        float theta_reflect = std::asin(
            1.0F / refract); //-------临界角--求low refractive index/high refractive index的arcsin

        if (theta_refract < theta_reflect)
        { // 小于临界角，要考虑吸收、反射、折射；大于临界角，要考虑TIR
            // 反射向量与玻璃面反面法线夹角小于全内反射角
            float eta_1 = refract;
            float k_1 = 1.0 - eta_1 * eta_1 * (1.0F - dot(-normal3, dir3) * dot(-normal3, dir3));
            //        if(k<0.0f)printf("k=%f \n",k);
            //        if(k_1<0.0f)printf("k1-%f,reflect=%f
            //        refract=%f\n",k_1,theta_reflect,theta_refract);
            ////如果结果有nan，测试这里
            float3 dir4 =
                normalize(eta_1 * dir3 - (eta_1 * dot(-normal3, dir3) + sqrt(k_1)) * (-normal3));
            //            float3 dir44=normalize(eta_1*dir333-(eta_1* dot(-normal3,dir333)+
            //            sqrt(k_1))*(-normal3));
            // printf("dir4 %f %f %f origin %f %f %f  dir44 %f %f %f origin1 %f %f
            // %f\n",dir4.x,dir4.y,dir4.z,origin1.x,origin1.y,origin1.z,dir44.x,dir44.y,dir44.z,origin13.x,origin13.y,origin13.z);

            // 折射向量
            //            printf("orgin %f %f %f now %f %f %f  normal %f %f %f now %f %f
            //            %f\n",origin.x,origin.y,origin.z,origin1.x,origin1.y,origin1.z,normal1.x,normal1.y,normal1.z,normal3.x,normal3.y,normal3.z);
            dir = dir4;
            origin = origin1;
            normal1 = normal3;
        }
        else
        { // 发生TIR
            // dir= make_float3(1.0,0.0,1.0);
            return;
        }
    }
    else
    {
        dir = normalize(reflect(-dir, normal1)); // 这里是直接进行反射，不考虑折射
    }
    // printf("dir %f %f %f \n",dir.x, dir.y, dir.z);
    // 计算经过定日镜的光线是否被遮挡----如果被遮挡，则不考虑与接收器相交的部分
    if (rectGridDDA::collision_(origin, dir, rectangleGrid, d_subHeliostat_vertexe,
                                heliostatArgument, no))
    { // ！！！！！！！！！！！！！！！！！！！！！！！！！！
        // printf("exit cuda --- collision-2\n");
        return;
    }
    // if(collision_in_rectgrid(origin, dir, rectangleGrid, d_subHeliostat_vertexe,
    // heliostatArgument,no)){
    //     return;
    // }

    // 计算光线是否能到达接收器
    //----尽管前面已经调整过定日镜的法向以将反射光对准接收器，但是由于太阳光到达定日镜是以光锥的形式，不是平行光，且定日镜表面有粗糙度会导致slopeerror，
    //----因此光线经过定日镜后，不一定能够与接收器相交
    float t;
    float u;
    float v; // t,u,v分别指的是什么
    if (!cylinderReceiver.GIntersect(origin, dir, t, u, v))
    {
        // printf("exit cuda --- miss rev\n");
        return;
    }
    //    printf("orign %f %f %f  dir %f %f %f  normal %f %f %f\n",origin.x, origin.y, origin.z,
    //    dir.x, dir.y, dir.z, normal1.x,normal1.y,normal1.z);
    // printf("t = %f, u = %f, v = %f\n",t, u, v);

    // printf("cuda factor = %f\n", factor);

    // 计算大气衰减
    float eta = 0.0;
    if (t <= 1000.0F)
        eta = 0.99331f - 0.0001176f * t + 1.97f * (1e-8f) * t * t;
    else
        eta = std::exp(-0.0001106f * t);
    float decay = 1.0F;
    float energy = fabsf(dot(dir, normal1)) * eta * factor * decay;
    //    printf("fabsf(dot(dir,normal1)) = %f, eta = %f, factor = %f \n", fabsf(dot(dir,normal1)),
    //    eta, factor );

    //    float energy=fabsf(dot(dir, normal1)) * eta* factor;

    // !!!!!!!!!这里对v的处理又是为什么？？？？？
    v = 1 - v;

    if (v < 0.25)
    {
        v += 0.75;
    }
    else
    {
        v -= 0.25;
    }

    // 把本次计算的光线的能量累加到接收器的像素上
    int nx = cylinderReceiver.resolution_.x;
    int ny = cylinderReceiver.resolution_.y;
    float* imagess = cylinderReceiver.getDeviceImage();
    int2 row_col = make_int2(u * ny, v * nx);
    int addresses = row_col.x * nx + row_col.y; // col_row.y + col_row.x * resolution.y
    // printf("u=%f, v=%f, energy=%f\n", u, v, energy);
    atomicAdd(&(imagess[addresses]), energy); // CUDA atomic

    //    printf("-------------helio_no = %d, micro_no = %d, energy = %f\n", helio_no, micro_no,
    //    energy);
}

void solar::cuda::CylinderReceiverRectangleGridRayTracing(SunrayArgument* sunrayArgument,
                                             CudaCylinderReceiver* cylinderReceiver,
                                             CudaRectangleGrid* rectGrid,
                                             HeliostatArgument* heliostatArgument,
                                             float3* d_subHeliostat_vertexe, float factor, bool res,
                                             solar::CloudLoader* cloud)
{
    printf("in func: CylinderReceiverRectangleGridRayTracing\n");
    int n_threads = 128;
    dim3 n_blocks;
    /*    printf("num of")

     *
     float mean_ =0.0022f;
     float err_ = 0.0015f;

     float *gaussian = new float[heliostat_nums];

     //Step 3: Generate theta and phi
     RandomGenerator::cpuGaussian(gaussian, mean_, err_, heliostat_nums);

     float *d_phi= nullptr;
     checkCudaErrors(cudaMalloc((void **) &d_phi,  sizeof(float) * heliostat_nums));

     RandomGenerator::gpuUniform(d_phi,heliostat_nums);


     float *d_gaussian= nullptr;
     checkCudaErrors(cudaMalloc((void **) &d_gaussian,  sizeof(float) * heliostat_nums));
     global_func::cpu2gpu(d_gaussian,gaussian,heliostat_nums);
 */
    // printf("%d %d \n",heliostat_nums,heliostatArgument.numberOfSubHeliostats);

    printf("about to enter func: setThreadBlocks\n");
    printf("heliostatArgument.helionums = %d\n", heliostatArgument->helionums);
    printf("heliostatArgument.micronums = %d\n", heliostatArgument->micronums);
    printf("sunrayArgument.numberOfLightsPerGroup = %d\n", sunrayArgument->numberOfLightsPerGroup);

    global_func::setThreadBlocks(
        n_blocks, n_threads,
        (long long)((long long)heliostatArgument->helionums * heliostatArgument->micronums *
                    (long long)sunrayArgument->numberOfLightsPerGroup),
        true);

    printf("num of rays= %lld \n", (long long)heliostatArgument->helionums *
                                       (long long)heliostatArgument->micronums *
                                       (long long)sunrayArgument->numberOfLightsPerGroup);

    printf("nBlocks = (%d,%d,%d), nThreads = %d\n", n_blocks.x, n_blocks.y, n_blocks.z, n_threads);
    printf("factor = %f\n", factor);
    map_raytracing<<<n_blocks, n_threads>>>(*sunrayArgument, *cylinderReceiver, *rectGrid,
                                            *heliostatArgument, d_subHeliostat_vertexe, factor, res,
                                            *cloud);
    cudaDeviceSynchronize();
    checkCudaErrors(cudaGetLastError());

    printf("----------cylinder_rev_rect_grid_raytracing_ends!!!-------\n");
}
